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FOREWORD 


The present report documents work performed under Contract 
NASl-10843 during the period October 1, 1979 to September 30, 
1980. 

The basic formulation and software for execution of the global 
function approach were developed under contract funding. In a 
separate but related effort, automatic procedures were developed 
for the selection of some of the parameters governing the anal- 
ysis strategy used in the program. This work was funded under 
the LMSC independent research program but is included in this 
report for completeness. 
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INTRODUCTION 


The structural response to a given environment is described by the differen- 
tial equations of motion of deformable bodies. Analytic solutions of such 
problems for a reasonably large class of structural configurations are not 
within the realm of the possible. Consequently, the mathematical problem 
is recast into a numerical problem for solution on the computer. 

The output from the computer consists of a sequence of numbers, in some way 
representing the functions satisfying equilibrium equations and boundary 
conditions. If the solution is represented by a linear combination of a set 
of "basis functions" then the components of the output vector consist of the 
coefficients in this linear combination. This is the case if we use the 
Galerkin or Rayleigh-Ritz procedures. If we use the finite difference or 
finite element procedures, the solution function is represented by its values 
at a number of discrete locations within the structure. Because these 
discretized methods are readily applied in a computer program for a general 
type of structure, they have been gaining popularity. This applies in 
particular to the finite element method. The finite element method may be 
considered as a Rayleigh-Rltz analysis in which the basis functions are 
localized. Confusion is avoided if the classical form of the Rayleigh- 
Ritz analysis is referred to as the "Global Function Approach". 

New technology in the space and energy fields has led to a growing demand for 
accurate analysis which at times cannot be met due to the limits set by 
available budget for computer time. In response to this need for more 
efficient numerical analysis, the possibilities have been explored of reduc- 
ing the number of freedoms in the system through a revival of the global 
function approach. Nagy (Reference 1) analyzed trusses using buckling modes 
as Ritz functions. 

This approach is straightforward if it can be assumed the deformation is 
inextensional, but is not directly applicable if the strain energy due to 
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stretching of the neutral axis (middle surface for shells) must be included. 

A more general way to automatically select a suitable set of basis vectors 
and a method to control the accuracy of the solution was first presented in 
Reference 2. In that case, the global function approach is used in connec- 
tion with a finite element model. The basis functions are represented by a 
set of basis vectors and the interpolating shape functions. Approximate 
solutions to the initial system are sought in the reduced space defined by 
all linear combinations of the basis vectors. We refer to the reduced space 
as the (infinite) set of trial vectors. Basis vectors can be obtained through 
solution of the initial (discrete) system and the accuracy of solutions in the 
reduced system can be assessed in the discrete system. Such procedures were 
further developed by Noor, References 3 and 4. In Reference 2, orthonormalized 
nonlinear solutions at different load levels are used as basis vectors while 
Noor proposed to use the so-called path derivatives. In both cases, the 
procedure involves a return to the discrete system (the finite element model) 
for evaluation of the error and automatic generation of new vectors when the 
size of the error suggests such action. 

In a finite-element formulation of the structural problem, we have 

MX + SX - F = 0 (1) 

where M is the mass matrix, S a nonlinear algebraic stiffness operator, 
and F the vector of external forces. The vector X represents the free- 
doms in the finite-element formulation; that is displacement and rotation 
components at the structural nodes. 

A global function formulation may be obtained by introduction into the 
finite element formulation of the substitution 


X = Tq 


( 2 ) 


where each column in the matrix T represents one of the basis vectors. 
A basis function is defined by the finite-element discretization, i.e.. 
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by the displacement and rotation components at structural nodes and by the 
local shape functions peculiar to the element. The components of the 
vector q are coefficients in a linear combination of basis vectors. These 
coefficients are the degrees of freedom in a reduced nonlinear algebraic 
equation system obtained through substitution of Eq. (2) into Eq. (1) and 
summation over the elements. The ith equation is of the form 


M.q. 

11 


+ I 


A. .q . 
iJ J 


lEB . q .q, + 




%kl^2 


^k 'll = ^1 


(i,j ,k,l = 1,1) 


(3) 


where M and F are generalized masses and forces corresponding to the 
1th basis function. The nonlinear terms derive from the stiffness operator. 


Solution Method 


An assemblage of computer programs was developed which is based on the use 
of global functions together with a finite element model. This assemblage 
includes the STAGSC-1 program (Reference 5) . The structural model is always 
defined by STAGSC-1 input data. In non-linear elastic analysis, the program 
user has options to define global functions (as input) or to obtain such 
functions through solution of the discrete system. Eigenmodes, buckling or 
vibration, or nonlinear solutions to the static equilibrium equations can 
be xncluded (users choice) . A special feature of the solution method is 
a problem-adaptive solution strategy in which automatic choice and 
continuous modification of certain strategy parameters allows for efficient 
analysis. 

The automatic feature is based on the ideas first proposed in Reference 2. 

In that case, the use of global functions represents a powerful way to deter- 
mine initial estimates for iterative solution of the discrete system. Initial 
estimates are obtained through Integration of a reduced displacement space 
spanned by the selected basis vectors. Successful operation of such software 
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requires the availability of adequate methods for specification of the 
basis vectors and a satisfactory step size selector sensing when a return 
to the discrete system for updating is desirable. In References 3 and 4, 

Noor uses the path derivatives as basis vectors. These are defined in terms 
of the coefficients in Eq. (3). Unfortunately, the number of distinct 
coefficients is very large and if many different elements are included in the 
structural model, severe disc storage problems will result. Therefore, the 

automatically selected basis vectors are defined in terms of nonlinear 
solutions to the discrete system at a sequence of different load levels. 

This is equivalent to use of numerically determined path derivatives. 

The disadvantage in this case is that solution accuracy may limit the 
number of solutions that profitably can be included as basis vectors. 

In Reference 3, Noor bases the return key on the change in a structural 
stiffness parameter. A change by ten percent in the value of this parameter 
prompts return to the discrete system. In Reference 4 he uses, as in 
Reference 2, the norm of the error in the equation system for the discrete 
model. In both cases, the analysis is very efficient. In Reference 4, the 
collapse analysis of an axially compressed prismatic shell (the "pear 
shaped cylinder") is presented. Return to the discrete system is dictated 
by an error norm (normalized with respect to the norm of the load vector) 
exceeding 0.05. The load steps i.e. the load increments between returns 
to the discrete system in that case are very large; indicating a potential 
for substantial savings in computer cost in non-linear elastic analysis. 


Extensive experimentation with automatic selection has indicated that both a 
stiffness parameter and the error norm may be useful for step size control. 
However, suitable values of the parameters governing return to the discrete 
system are not only case dependent but in a given case may may also vary 
considerably with the load level. It appears that the potential for savings 
in computer time by use of global functions can only be realized if a problem- 
adaptive computational strategy is available. Then, the return key is adjusted 
in response to the characteristics of the problem so that an efficient 
analysis can be obtained in a variety of cases without preceding 
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experimentation with the step size selection. The procedure involves a number 
of strategy parameters. A set of default values for these parameters has been 
selected to be used when the analyst lacks special knowledge of the behavior 
of his system and therefore declines to make a different choice. 

Default Strategy 

The initial basis vectors are obtained through solution of STAGSC-1. This 

program does not contain procedures for automatic choice of the initial load 

step. The user must define the ititial load and the initial step in 

STAGSC-1. The user also defines the initial and the maximum numbers of 

basis vectors, N, and N . Default values are N. = 4, N =6. 

1 m 1 m 

On each return to the discrete system, a nonlinear solution of this discrete 
system is obtained and included in the data base. If the number of basis 
vectors in the data base exceeds N^, the program gives preference to those 
corresponding to higher load levels when the basis vectors are selected. 

A check on linear dependence among the basis vectors is performed and vectors 
that are not sufficiently distinct are discarded. The number of basis vectors 
therefore can be less than and remain less than N^. During the computa- 
tions, the program attempts to set the return key so that solution of the 
discrete system will require approximately five iterations. The stiffness 
parameter included in the strategy is represented by the diagonal elements 
in the factored matrix corresponding to the reduced system. The following 
notations are used: 

e = error norm = 1 1 SU | | / ] ] f j ] where 6U is the first variation 
of the total potential energy (i.e., the residuals) and f 
is the vector of applied forces (including reactions) . 

A = vector of diagonal elements in the factored matrix 
associated with the reduced system. 

A. = ratio between present value of the elements of A and 

corresponding values at last return to the discrete system 
(initial values for the first step) . 

N = number of iterations for convergence at last return to the 
discrete system. 

-k(N-5) 

f = 10 ^ ' where k is an input constant. 
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Solution o£ the discrete system and updating of the set of basis vectors by 
inclusion of the current solution is dictated by any of the following events: 


(1) e > 

(2) for some i 

(3) A. < 6. for some i 

(4) AP is reached (max load step) 

where c , 6„, 6^ are input parameters, 
q U L 

Whenever convergence occurs on the return to the discrete system, the 
adjustment depends on which criterion prompted the return. 


If e > e , then e f e 

q q q 

A. > A,,, then 6 -»• 1 + f (<S„ - 1) 

1 U’ U U ^5^ 

A. < 5 , then 6. ->■ 1/[1 + f (l/'^-r ~ 1)] 

1 L L t 

AP reached, then AP -»■ f AP 

If divergence occurs on return to the discrete system, f is set to 0.5 and 
all four return keys are accordingly adjusted. 

Initial efforts established that efficient analysis would be achieved in a 
variety of cases with the choice: 

e = 0.2, =4,6= 0.3, and k = 0.08 (6) 

q U 1j 

strategy should be considered as a first cut only • Additional iinprove” 
ments are certainly possible. For example, a good initial value of can 
probably be surmised from the relation between error norm and convergence 
rate in the first series of solutions in STAGSC-1. Also it may be better 
to adjust all the strategy parameters on any return to the discrete system, 
while a more efficient automated strategy may be forthcoming, the present 
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version was evaluated through comparison to solution with STAGSC-1 in a 
study of five structural configurations with significant nonlinearity. 
The automated strategy used was in all cases based on the default values 
for the parameters. 


Benchmark Cases 

The automated global function procedure is primarily intended for use in the 
analysis of cases requiring large amounts of computer time. In such cases, 
the time spent on integration in the reduced system would be insignificant. 
Further reductions of the computer time in auxiliary programs is possible, 
that is, in the routine used to define and to integrate the reduced system. 

In terms of required computer time, the five benchmark cases range somewhat 
below the class of problems for which computer cost is a serious issue. A 
straight comparison of computer run times required by the global function 
analysis and the standards STAGSC-1 may not be conclusive. Therefore, the 
total number of iterations and the number of reformulations and factorings 
of the second variation are also recorded in each case. In all the cases, 
the analysis was carried well into the non-linear range but never beyond 
possible limit points since STAGSC-1 presently does not include efficient 
methods to handle such cases. 

Case 1; "Pear-Shaped Cyliner" in Compression 

The first test case is the Pear-Shaped Cylinder considered in Reference 7. 
Dimensions and material data are shown in Figure 1. The cylinder is subjected 
to uniform axial shortening. The maximum load factor was set to 8.27 correspond- 
ing to an axial shortening of 00.00004201 m. which is within a couple of percent 
of the critical value. The corresponding axial load is 12793 N. This value is 
well above the critical load reported in Reference 6, possibly caused by a 
tendency of the finite element configuration used to "lock" as the rotations 
become relatively large. Whatever the reason, this discrepancy has no impact 
on the solution procedures. 


8 



Fig. 1 Deformation of Pear-Shaped Cylinder 
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Due to the presence of two syimnetry planes, only one-quarter of the shell 
need be considered. A finite difference grid is defined over this model with 
5 uniformly spaced gridlines axially and 43 approximately uniformly spaced 
gridlines circumferentially. 

While the analysis with STAGSC-1 required a total of 226 iterations and 13 
refactorings, the corresponding numbers with the automated global function 
analysis are 25 iterations and 4 refactorings. The total run time (CP-time 
CDC 175,NOS-BE) is 605 seconds with STAGSC-1 and 240 seconds with the global 
analysis. However, with 1300 degrees of freedom and an average bandwlth of 
40 this is a relatively inexpensive case. The time spent on formulation 
and solution of the reduced system is still significant, approximately 
half the total run time. The indication then is that for a larger case the 
use of global functions for extrapolation may reduce the run time by a factor 
between 2.5 and 5.0. 

The load level at which STAGSC-1 solutions were obtained are indicated in 
Figure 1. 

The first three solutions represent results from the initial STAGSC-1 analysis 
The first return was dictated by excessive change in the stiffness parameter, 
the second by the size of the error norm. In each of these instances, the 
number of iterations in the discrete systems is five and no parameter adjust- 
ment was required. On the last return (maximum step size) the number of 
iterations is six. 

Case 2; Bending of a Long Cylinder 

As a second test case, a long cylinder subject to a constant bending moment 
was selected. The cylinder has a radius of 0.127 m. and a length of 3.048 m. 
as shown in Figure 2. Due to the symmetry conditions, only one-quarter of the 
cylinder needs to be included in the structural model. Over this model a 
uniform 11 x 11 finite element grid is applied. This corresponds to a total 
of 939 degrees of freedom and an average bandwidth of 82. The bending moment 
is applied at one end and symmetry conditions are enforced at the other. 
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Fig. 2 Bending of Long Cylinder 
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A hefty end ring is applied at the loaded end, preventing significant cross- 
section warping. 

Since convergence of the buckling load is from above, the shortwave local 
buckling observed in Reference 7 is suppressed by use of a coarse net and a 
pure Brazier effect is displayed. For an infinitely long cylinder. Brazier, 
Reference 8, predicts collapse under a bending moment of 

w_2/? Eirah^ ,, „ 

M ^ = 14,620 N-m (7) 

/ 1 -v^ 

The cross-section at the point of collapse is then flattened by 0.05588 m. 

(0.22 times the diameter). 

The results of the STAGSC-1 shown in Figure 2 appear to be in relatively 
good agreement with those by Brazier. The analysis was interrupted when 
the bending moment reached 13445 N’m. At that point, convergence is very 
difficult and the load step must be very small. With more accurate extrapola- 
tion, the global function analysis was carried somewhat further 13705 N*m. 

The points at which STAGSC-1 solutions were obtained (including the initial 
solutions and subsequent returns for updating) are indicated in Figure 2 
together with two points at which convergence failed. The total run time 
with the global function analysis was 335 sec CP as compared to 1030 with 
STAGSC-1. With the global analysis, the total number of iterations is 50 
and the number of refactorings 20. Corresponding numbers with STAGSC-1 are 
450 and 34. In this case, about 75 percent of the global analysis run time 
is spent in STAGSC-1. The indication for a large case with this general 
behavior is a reduction in computer cost by use of global function extrapola- 
tion of 3 to 4. 

Each return to the discrete system was caused by excessive error in the first 
variation of the discrete system. Figure 3 shows how this error varies with 
the applied load. At each point of update, the error returns to zero. The 
default value = 0.2, is for this case somewhat too small; convergence 
is obtained in two iterations and the value of e is eventually increased 
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e 

Cq 

1.0 — ^ 

O NUMBER OF ITERATIONS 

-f INDICATES FAILURE TO CONVERGE IN DISCRETE SYSTEM 

e,eq 




to 0.6. Convergence becomes more difficult as the critical load is approached 
and c IS reduced by a factor of four after divergence in two consecutive 

q 

efforts . 

Case 3: Spherical Shell Subjected to Point Force 

The case of a spherical cap clamped at the edge and subjected to an inward 
directed central point force at the apex is considered as a third case. The 
shell geometry shown in Figure 4 is identical to that considered by Fitch in 
Reference 10. The bifurcation buckling analysis with nonlinear prestress in 
Reference 10 indicated buckling into a mode with four circumferential waves 
when the displacement under the point load reaches 0.046736 m. 

The case considered here is the nonlinear behavior of the cap when in addition 
to the central point load two small forces are applied such that deformation 
in a four-wave pattern is triggered. These forces are held constant and the 
midpoint displacement is gradually Increased. A 90-degree sector of the cap 
is analyzed with symmetry conditions applied along the meridians. The finite 
difference grid is uniform in the circumferential direction and varied in the 
meridional direction with the grid spacing tightened around the apex and at 
the clamped edge. The grid includes 16 grid lines in each direction. This 
system contains 1966 degrees of freedom and the average bandwith is 127. 

The STAGSC-1 works quite well up to a midpoint displacement of about 0.0254 m. 
After that convergence is more difficult. The computations were interrupted 
when the displacement reached 0.040132 m. This does not correspond to a 
maximum in the load displacement curve but even with a step size of 0.000254 m. 
(deformation) it is at that point necessary to refactor on each load step. 

The results including the growth of a four-wave buckling pattern are shown 
in Figure 4. 

With the global function analysis the maximum load factor, 2.0, was easily 
reached. Load levels at which STAGSC-1 solutions were obtained are shown 
in Figure 4. The error bound, c = 0.2, is much too severe in this case 


14 




Fig. 4 Deformation of Spherical Cap 
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the AUTORITZ analysis Is not efficient. However, the return key is adjusted 
(to = 10.5) and a suitable step sxze is eventually chosen. A displacement 
of 0.040132 m. was reached with a total of 52 iterations and 20 reformulations 
with factoring of the second variation. The corresponding values with STAGSC-l 
are 196 iterations and 20 refactorings. The run time is 1898 sec with STAGSC-l 
and 971 sec with AUTORITZ of which about 70 percent is spent in the STAGSC-l 
program. The saving in computer time for a case of this type then is 2.0 to 
2.7 and it will be considerably more in favor of AUTORITZ if the analysis is 
carried further. 


Case 4; Panel with Initial Imperfections 


The fourth case is a cylindrical panel subjected to axial compression in the 
form of uniform end-shortening. The properties of the panel are shown in 
Figure 5. All the four edges are simply supported. The panel is free from 
initial stresses but deviates from the true geometric shape. The initial 
imperfection is represented by a lateral displacement (in meters) of the form 


= 0.00127 sin ^ sin (12 y) + 0.000508 sin ^ sin (24 y) 
+ 0.000508 sin sin (12 y) + 0.000508 sin sin (24 y) 

J-i J-i 

+ 0.000508 sin sin (12 y) + 0.000254 sin ^ sin (24 y) 


( 8 ) 


+ 0.0000508 sin ^ sin (36 y) 

Ij 


A uniform 19 x 19 finite element grid was used resulting in a system with 
3044 degrees of freedom and an average bandwidth of 151. The lateral 
displacement at the midpoint and the axial shortening are shown in Figure 5 
as functions of the total axial load. With STAGSC-l the maximum load, 
3,247,200 N. , corresponding to an axial shortening of 0.01143 m. , is reached 
after a total of 106 iterations and 9 refactorings. The total run time is 
1096 sec. The corresponding values with AUTORITZ are 23 iterations, 

6 refactorings, and a total run time of 786 sec. The three initial load steps 
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Fig. 5 Deformation of Imperfect Panel 


and the returns to the discrete system for updating are indicated in the 
figure. The first three returns are governed by the size of the maximum load 
step. This is the case in which the global analysis compares least favorably 
with a straight STAGSC-1 analysis. The reason seems to be that the STAGSC-1 
strategy works unusually well in this case. The saving in computer time by 
use of global functions in this case is only by a factor of 1.4. With 
further improvement of the efficiency of the automated strategy, a saving 
by a factor of two seems possible. 

Case 5; Cylinder with Cutouts 

The last case considered was a cylinder with two diametrically opposite 
rectangular cutouts. The geometry of the shell together with some results 
of the analysis are shown in Figure 6, Here, u^ represents the uniform end 
shortening and W the lateral displacement midways on the edge of the cut- 
out. As in Case 1 (pear shaped cylinder), the maximum load reached exceeds 
the previously established in analysis as well as in experiments (Reference 
12) . This tendency of the element to lock with relatively large rotations 
occurs despite the fact that a very fine grid has been used as indicated in 
Figure 7. The model has 7055 degrees of freedom and the average bandwidth 
is 205. 

In the analysis, the uniform end shortening was gradually increased to 
0.00009144 m. which is just below the limit point for the model and corresponds 
to a total force of 12900 N. With STAGSC-1, this load level was reached with 
a total of 145 Iterations and 12 refactorings. By use of the automated 
strategy, it was only necessary to formulate and factor the second variation 
6 times and the total number of iterations was 16. The total run time with 
STAGSC-1 was 5267 sec and the global analysis 1856 sec CP time; hence, the 
use of global functions for extrapolations leads to a considerable saving - 
a factor 2.8 between the run times. 
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Fig. 6 Deformation of Cylinder with Cutout 
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Fig. 7 Finite Element Grid for Cylinder With Cutout 


20 



CONCLUSIONS 


The present paper discusses the possibilities to improve the efficiency in 
nonlinear structural analysis by use of global displacement functions. 

A procedure is presented for accuracy control and automatic generation of the 
global functions. In essence, the procedure is a way to predict initial 
estimates for solution of large nonlinear systems with a continuation method. 
The indication from a few benchmark cases is that an improvement by a factor 
of three to five is possible in most cases. The sample cases indicate the 
usefulness of the procedure in solution of nonlinear structural shell problems 
by the finite element method. However, the basic idea of extrapolation 
through integration in a reduced solution space may be useful in other applica- 
tions leading to large and strongly nonlinear equation systems. 
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